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Abstract 


The inverse scattering problem for an acoustic medium is formulated by 
using the variable background Born approximation. A constant density acoustic 
medium is probed by a wide-band point source, and the scattered field is 
observed along a curved receiver array located outside the region where the 
medium velocity is different from the assumed background velocity function. 

The solution that we propose relies on the introduction of a backpropagated 
field. This field is obtained by using a finite-difference scheme backwards 
in time to backpropagate into the medium the scattered field observed along 
the receiver array. The backpropagated field is imaged at the source travel 
times, giving an image of the same type as obtained by reverse-time 
finite-difference migration techniques. The gradient of this image is then 
taken along rays linking the source to points in the medium, and after 
scaling, this gives the reconstructed potential. To relate the reconstructed 
potential to the true scattering potential, we use high frequency asymptotics 
and an additional approximation introduced by Beylkin [1]. These 
approximations reduce the validity of our reconstruction procedure to the high 
wavenumber region. With these approximations, it is shown that at a given 
point, the reconstructed potential corresponds to the convolution of the true 
potential with a weighting function obtained by partially reconstructing an 
impulse from its projections inside a cone. The angular range of this cone is 
totally determined by the geometry of the receiver array, and by the relative 
location of the source with respect to the point that we consider. In the 
special case when the receiver array surrounds the domain where the scattering 
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potential is located, we find that within the Born approximation, the 
reconstructed potential recovers exactly the high wavenumber part of the 
Fourier transform of the true potential. It is expected that for a wide class 
of problems, the reconstruction technique described in this paper will be 
computationally more efficient that the generalized Radon transform (GRT) 
inversion method proposed by Beylkin, Miller and Oristaglio [1] - [3]. 
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1. Introduction 

The multidimensional inverse scattering problem for acoustic media has a 
wide range of applications in areas such as exploration geophysics [4], 
ultrasonic imaging [5], and nondestructive testing, among others. In this 
problem, the medium of interest is probed by sources located outside the 
medium and the scattered field is recorded at various locations. The 
objective is to reconstruct the velocity and density of the medium as 
functions of position. In this paper, we will assume for simplicity that the 
density of the medium is constant. Since the relation between the 
observations and the velocity function is nonlinear, the solution of the 
inverse problem must also be nonlinear. For one-dimensional (1-D) media, i.e. 
for media whose velocity function varies only with one space dimension, a 
number of exact inverse scattering procedures which rely either on integral or 
on differential equations have been proposed over the years. A discussion of 
these methods can be found in [6]. However for multidimensional media, exact 
inverse scattering methods are still at an early stage of development, and the 
methods which have been developed until now require experiment geometries 
which are rather unrealistic from a practical point of view. 

This has motivated researchers to develop approximate direct inversion 
methods, which solve exactly a linearized form of the multidimensional inverse 
scattering problem. The Born and Rytov approximations [7] are the most 
commonly employed of these linearization techniques. Whether one should use 
one of these approximations instead of the other depends on the nature of the 
scatterers and on the experiment geometry. In this paper, we shall consider 
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the Born approximation. In this approximation, the object profile that we 
want to reconstruct is viewed as a small perturbation about an assumed 
background velocity model, and the scattered field is expressed linearly in 
terms of this perturbation. Physically, the Born approximation takes into 
account only the singly scattered waves; multiple scattered waves due to the 
velocity perturbations are neglected. Note however that the multiples due to 
the background model are included in the scattered field. 

During the past few years, a number of solutions of the multidimensional 
Born inversion problem have been proposed for various observation geometries 
and background velocity models. The majority of these solutions assume a 
homogeneous background model. The zero-offset reflection geometry consisting 
of coincident sources and receivers was considered by Cohen and Bleistein [8] 
for a line aperture in two dimensions, and by Norton and Linzer [5] for plane, 
cylindrical and spherical apertures in three dimensions. More recently, 
Fawcett [9] formulated the zero-offset Born inversion problem as a tomographic 
problem, where the objective is to reconstruct a function from its projections 
along circles or spheres. The inversion method proposed by Cohen and 
Bleistein was also extended to the case of a stratified (1-D) background model 
in [10] and [11]. Raz [12] and Clayton and Stolt [13] considered the same 
experiment geometry as Cohen and Bleistein, but with unstacked data. In this 
geometry, some source and receiver arrays are located on the surface of the 
earth, and for each source the scattered field is recorded at all receivers 
rather than only at the coincident receiver. They showed that both the 
density and bulk modulus of the acoustic medium can be recovered with this 
experiment. Another interesting feature of Clayton and Stolt’s paper is that 
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it assumes a variable background model, which is then taken into account by 
extrapolating the observed scattered field downwards into the medium of 
interest. 

The solution of the multidimensional Born inversion problem which was 
proposed by Esmersoy and his colleagues [14]-[16] relies on a similar 
backpropagation principle. However, they considered a different scattering 
experiment, where the medium is probed by a single wide-band point or 
plane-wave source, and where the scattered field is measured along a curved 
receiver array located outside the region where the medium velocities differ 
from the background model. In this approach, the scattering potential is 
reconstructed by propagating the observed scattered field backwards in time 
with a finite-difference scheme, imaging this field either at zero time [15] 
or at the source travel times [14], [16], and then filtering the resulting 
image. Depending on whether the receiver array surrounds the scattering 
object or not, the reconstructed potential is an exact or approximate solution 
of the linearized inverse scattering problem. However, an important 
limitation of the results described in [14]-[16] is that they were restricted 
to the homogeneous background case. 

The first complete solution of the variable background Born inversion 
problem was presented by Beylkin, Miller and Oristaglio [1] - [3], who 
formulated this problem as a generalized Radon transform (GRT) inversion 
problem, where the objective is to reconstruct a function from its projections 
along a set of curves whose geometry depends on the experiment and on the 
background model. The solution obtained by Beylkin and his colleagues is 
expressed in terms of a weighted backprojection operator slimming the 
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contributions of the curves passing through a given point. This solution is 
quite general since it applies both to the case when the medium is probed by a 
single point source, and to the zero-offset and finite-offset geometries. 

Note however that the inverse GRT relies on high frequency asymptotics, and on 
an additional approximation which together have the effect that the 
reconstructed potential recovers only the rapid space variations of the true 
potential, or equivalently the high wavenumber part of its Fourier transform. 
Thus, the inverse GRT solves the variable background Born inversion problem 
only in a limited sense. Another limitation of the.inverse GRT is that to 
obtain the weights appearing in the backprojection operator, one must use a 
ray tracing algorithm for every point along the receiver array where 
measurements are taken. The amount of computations required by this method is 
therefore quite large. 

The objective of this paper is to extend the backpropagated field 
approach of [14]-[16] to the variable background case for the experiment where 
the medium is probed by a single wide-band point source. The computational 
procedure that we use to obtain the reconstructed potential is quite different 
from the inverse GRT, but the domain of validity of our reconstruction method 
is the same as that of the inverse GRT. In other words, the reconstructed 
potential recovers only the rapid space variations of the true potential, or 
equivalently the high wavenumber part of its Fourier transform. The first 
step in our approach is to filter the observed time traces and to use them as 
source wavelets for doublet sources located along the receiver array. For 
this source configuration, a finite difference scheme is used backwards in 
time to compute the backpropagated field inside the medium. This field is 
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then imaged at the source travel times, and the reconstructed potential is 
obtained by taking the gradient of this image along rays linking the source to 
points in the medium, and by scaling the resulting expression. To relate the 
reconstructed potential to the true potential, we use high frequency 
asymptotics as well as an additional approximation introduced by [1], These 
approximations reduce the domain of validity of our inversion method to the 
high wavenumber region. In this context, it is shown that at a fixed point, 
the reconstructed potential is the convolution of the true potential with a 
weighting function obtained by partially reconstructing an impulse from its 
projection inside a cone. The angular range of this cone is totally 
determined by the geometry of the receiver array and by the location of the 
source with respect to the point that we consider. A consequence of this 
representation is that in the special case when the receiver array provides a 
total coverage of the scattering region, and within the domain of the validity 
of the Born approximation, our reconstruction procedure recovers exactly the 
high wavenumber part of the Fourier transform of the scattering potential. 

The step of our reconstruction method which is the most demanding from a 
computational point of view is the computation of the backpropagated field 
with a finite-difference scheme. We expect that for a significant number of 
problems, this computation will be easier to perform than the ray tracing 
scheme which is required for every receiver by the inverse GRT. 

This paper is organized as follows. In section 2, the velocity inversion 
problem is formulated within the Born approximation. The definition and 
implementation of the backpropagated field are discussed in section 3. Our 
reconstruction method is described in section 4, and a representation theorem 
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is obtained for the weighting function relating the reconstructed potential to 
the true potential. This representation relies on approximations which reduce 
the domain of validity of our reconstruction procedure to the high wavenumber 
region. The representation that we obtain shows that at a given point, the 
weighting function is the partial reconstruction of an impulse from its 
projections over a cone. The construction of this cone is examined in section 
5, and is shown to depend only on the experiment geometry. Section 6 contains 
some conclusions and some thoughts about further extensions of our results. 
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2. Problem Description 

Consider the scattering experiment described in Fig. 1. A constant 
density 2-D acoustic medium is probed by an impulsive point source located at 
r 7 in the x-y plane. Since this 2-D medium is in fact a 3-D medium whose 
velocity function c(x) does not vary with the third dimension z, the 2-D point 
source is in fact implemented by a line source parallel to the z axis. The 
scattered field is observed along a curved receiver array T, and in the 
following it is assumed that T is a smooth curve parameterized by the arc 
length s, sel. 

Note that since we assume that the 2-D medium is probed by a line source, 
the experiment geometry described above is not totally realistic. In 
practice, it is much cheaper to use a point source to probe the medium. In 
this case the problem becomes a 2 1/2-D problem in the sense that although the 
velocity function c(x) varies only with two space dimensions, the waves 
propagate in three dimensions. The geometrical spreading associated to the 
three-dimensional propagation of the waves needs to be taken explicitly into 
account in the inversion problem, and a detailed analysis showing how this can 
be done can be found in [17], and [14], section 2.6. However, to simplify our 
analysis, it will be assumed below that the medium is probed by a line source. 

Then, the Fourier transform P(x,u) of the pressure field satisfies 

[v 2 +k 2 n 2 (x)] P(x,w) = -S(x-rj) (2.1) 

where k = w/c is the wavenumber, and n(x) = c/c(x) is the refraction index of 
the medium measured with respect to some constant velocity c. Throughout this 
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paper, it will be assumed that n(x) does not deviate significantly from a 
background index n^(x), so that 

n 2 (x) = n o(*) + u (*) • C 2 - 2 ) 

where the scattering potential U(x) is small. In addition, we assume that 

Uq(•) is a smooth function and that U(•} has a bounded support V, which is 

completely located on the same side of T, as shown in Fig. 1. Substituting 

(2.2) inside (2.1), we can rewrite equation (2.1) as 

[v^ + k^n^fx)] P(x,cj) = - 5(x-ij) - k^U(x) P(x,w) (2.3) 

2 2 2 

where the operator Dq = v + k n^(x) appearing on the left-hand side of (2.3) 
is the background Helmoltz operator, and where the second term on the 
right-hand side can be viewed as a perturbation. The solution of (2.3) is 
given by 


P(x,u) = P 0 (x,w) 


+ k 


| dx'U(x' )P(x' ,w)Gq(x,x' ,gj) , (2.4) 


where the incident field Pq(x,u) satisfies the unperturbed equation 

Vo^") = C 2 - 5 ) 


and where Gq(x,x' ,«) is the Green’s function associated to Dq, i.e. 


DqGq(x,x',«) = -6 (x-x') 


( 2 . 6 ) 


By comparing (2.5) and (2.6) we see that for the experiment geometry 
considered here, the incident field Pq(x,w) = G^x.^.w). 

The main feature of equation (2.4) is that it is exact , i.e. no 
approximations are involved up to this point. This equation is known as the 
Lippmann-Schwinger equation [18], and it puts in evidence the nonlinear 
relation existing between the potential U(x) and the pressure field P(x,«). 

In this paper, we will linearize this equation by using the Born 
approximation , whereby the total field P(x,w) is approximated by the incident 
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field P Q (x, u) = Gq(x, tj, w) inside the integral appearing in equation (2.4). 
In this case, the scattered field P s (x,«) = P(x,u) - P Q (x,w) can be expressed 
as 


P g (x,(o) = k 2 J dx'U(x' )Gq(x,x‘ ,w)Gq(x' ,rj,<o) . (2.7) 

Note that the Born approximation is valid only if the scattering potential 
U(*) is small both in magnitude and spatial extent ([19], Chapter 9). Then, 
the relation between P g (*,w) and U(*) becomes 1inear . and the inverse problem 
that we shall consider in this paper consists in reconstructing U(x) for x a V 
from the observed scattered field P g (£,w) where £ = £(s) is located along the 
array T. Note that the receiver array T may not provide a total coverage of 
the domain V, so that in general U(x) will only be partially reconstructed 
from the given observations. 

Throughout this paper, we shall use the geometrical optics approximation 

e iir/4 

~ 1/2 a ( 5 . 5 , ) e lM5> - ^ ( 2 . 8 ) 

k 1 


for points x e V, and for x' = tj or x' a F, and where for negative values of 
1/2 1/2 

k, k = i(-k) . This approximation is a high-frequency approximation. 

Equivalently, it corresponds to assuming that the distance between the domain 
V where the scattering potential U(») is located and the source and receivers 
is large when compared to the wavelengths used to probe the medium, so that 

k|x-x’| » 1 (2.9) 

for x e V and x’ = tj or x’e T. In (2.8) ^(x.x’) is the phase or travel time 
function and satisfies the eikonal equation 


l v x f > (x*x') I = n Q (x), 


(2.10) 
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and the amplitude a(x,x‘) obeys the transport equation 


av <p + 2 v a*v <p = 0 

X X X 


( 2 . 11 ) 


along the ray linking points x and x'. For the special case when the 
background refraction index is constant, i.e. when n^x) = n^, we have 
GqCx.x',cj) = iH^^knQlx-x' |)/4 where denotes the Hankel function of 

order zero and type one, and the approximation (2.8) reduces to 

1 


iir/4 


““ITT — —i— 77J/2 

k 2(27 tt1qI x-x I) 


ikn^lx-x'| 


( 2 . 12 ) 


Then, substituting the approximation (2.8) inside (2.7), and setting 
x = £, we obtain 


where 


w) = f dx'U(x' )A(x‘ 

J v 

(2.13) 

A(x,f) = a(x,f)a(x,] 3 ) 

(2.14a) 

$(x,£) = <p(x,£) + <#>(x, 2 ) . 

(2.14b) 

F(£.k) = P s (£,co)/ik 

(2.15) 


Denoting 


and letting f(£,r) be the inverse Fourier transform of F(£,k) with respect to 
k, we find that 


f(f.r) = - [ P s (£,T)dr 

J° s 

f dx 1 U(x 1 )A(x',£)6(r-<J>(x\£)) 
Jv 


(2.16) 


where the time function P g (£.t) is the scattered field observed at £. 

The identity (2.16) indicates that f(£,r) can be viewed as a weighted 
projection of the potential U(•) along the curve $(x,£) = r, and that f(£,r) 
is obtained by integrating the scattered field p g (£,«). Thus, the inverse 
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scattering problem that we consider here can be formulated as a generalized 
Radon transform (GRT) inversion problem. This was the point of view adopted 
by Beylkin, Miller and Oristaglio [l]-[3], who were then able to obtain a 
backprojection operator to reconstruct the potential U(») partially from its 
projections f(£,r) where £ e T and 0 < r < °°. However, one disadvantage of 
the GRT inversion method is that to implement the backprojection operator, it 
requires the computation of the amplitude a(x,£) and phase <f>(x,£) for every 
point x in the medium and every receiver f e T, as well as the computation of 
a(x,Tj) and <#>(x, rj) for all points x, where the position rj of the source is 
fixed. This amount of computations is very large, and the inverse GRT is 
therefore very costly to implement. 

The objective of this paper is to obtain an alternative reconstruction 
method, whose computational requirements will be smaller for a significant 
class of problems. 

3. The Backpropagated Field 

The inversion method that we propose relies on the concept of 
backpropagated field, which was first introduced in the context of holographic 
imaging by Porter [20]. This idea was subsequently used by Bojarski [21] and 
Esmersoy [14] to study inverse source and inverse scattering problems, and it 
was applied to the solution of the constant background inversion problem for 
an impulsive point source and for a plane wave source in [15] and [16], 
respectively. The migration methods [22] - [24] which are currently used in 
exploration geophysics also rely on the concept of extrapolated field to image 
the main reflectors contained in a scattering medium. By migration we mean 
here a technique which is used to image the discontinuities of the scattering 
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potential U(x). Migration differs from inversion by the fact that in 
migration the objective is only to detect discontinuities of U(x), whereas 
inversion methods seek to obtain some precise quantitative information, 
perhaps only partial, about the values of the function U(x) or of its Fourier 
transform. 

For the experiment geometry considered here, the backpropagated field is 
defined as 

P e (x,co) = J ds W(«)P s (f.o,) V£G 0 *(x,£, W ).r(£) (3.1) 

where W(w) is a filter to be determined later, and where u(£) is the unit 
vector perpendicular to T at point £, oriented in the outwards direction, as 
shown in Fig. 2. Note that in integral (3.1), the receiver position £ = £(s) 
is a function of the arc length s, but for simplicity this dependence will not 
be indicated explicitly in the mathematical expressions that we shall consider 
below, except at places where our analysis will need to take this dependence 
into account. Here, v^Gq(x,£,u)*u(£) denotes the Green’s function of a 

doublet, i.e. it can be implemented by two impulsive sources of opposite signs 
located perpendicularly across the curve T at a very small distance from each 
other, as shown in Fig. 2. The fact that we select the complex conjugate of 
v^Gq(x,£,o)) indicates that the waves propagate anticausallv (backwards in 

time), since we want to reconstruct the pressure field inside the medium at 
earlier times. 

The expression (3.1) for the backpropagated field P e (x,u) has the 
following physical interpretation: P e (x,w) is the acoustic field which is 
obtained by replacing the receivers along the curve f by doublet sources, and 



16 


by exciting the doublet located at point £ e T with the source wavelet 

W(w) P g (£ .<*>)• This source wavelet is obtained by passing the scattered field 

observed at £ through the filter W(w). From a practical point of view, 

P e (x,w) can be computed in a number of ways. The method that we propose 
consists in observing from (3.1) that P e (x,to) satisfies 

D 0 P e (x,<o) = 0 (3.2) 

for x € T. Then if we specify an appropriate boundary condition for this 
equation, the backpropagated field P e (x,t), where P e (x,t) denotes the inverse 
Fourier transform of P e (x,w), can be computed backwards in time by using a 
finite-difference scheme of the type described in [25]. Note that except for 
the introduction of the filter W(w), the procedure described above for 
computing p e (x,t) is identical to the reverse-time finite difference migration 
method for unstacked data which was proposed in [26], where the boundary 
condition 

P e (£.t) = P s (£.t) (3.3) 

for [ e T and t > 0 was selected. 

The only element that has been left unspecified above is the choice of 
boundary conditions for p e (x,t). Strictly speaking, to obtain an exact 
boundary condition for P e (x,u) one would need to use the expression (3.1) to 
specify the value of P e (x,t) over a boundary F located close to the array F. 
However, it is more convenient to note that when equation (3.1) for P e (x,w) is 
equated to expression (3.11) below, and when T surrounds the domain V, it was 
shown by Esmersoy in [14], pp. 99-105 that 

p„(£.t> = pJ(E.t) 


(3.4a) 
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for £ £ T and 


t > max (<p(x, tj) - <p(x,£)) , (3.4b) 

xeV 

f 

where p g (£,t) denotes the inverse Fourier transform of the filtered scattered 
field W(w)P s (£,w). 

The identity (3.4a) specifies a boundary condition for the extrapolated 
field P e (x, t) over the time window (3.4b). Since the inversion method that we 
propose in section 4 requires the knowledge of P e (x, t) only at t = <#>(*.I?). 
where <#>(x.t]) is the travel time from the source 3 to point x e V, the time 
window (3.4b) is in general sufficient to compute the values of P e (x, t) that 
we need. In addition, the boundary condition (3.4a) is so simple that even 
when T does not surround the domain V, it may be worth using it to compute an 
approximation to P e (x, t). This scheme has given good results in [14], [16]. 
However, as indicated above, a more rigorous method would need to start from 
the definition (3.1) of P e (x,w). 

A representation of the backpropagated field which will be useful in 
subsequent derivations can be obtained by noting that within the geometrical 
optics approximation 

iw/4 ( 

v^G 0 (x,g, W ) = t - (v £ a(x,£) + ika(x,£) Y>(x,£))e lk ^’£ j (3.5) 

kl/2 

Then, if we make the additional approximation 

|v^a(x,g)| « |ka(x,£)| ( 3 . 6 ) 

which is of the same nature as the geometrical optics approximation, we find 


that 


7 £G 0 (x,£,(o) * -k 1 / 2 e ' i7r/4 a(x,£) v ^(x,f)e ik<p ^ 1 ^ ) 


(3.7) 
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A A 

The vector v(x,£) = Vg</>(x,£) is tangent to the ray linking x and £, as shown 
in Fig. 2, and 

|v(x,£)| = n 0 (£). (3.8) 

so that 

V £ G 0 (*’£’ W )* U (£) - “ k 1/2 e lir/4 n 0 (f)cos/3(x,g)a(x,f)e lk<p ^-'^^ 

= ikn 0 (£)cosP(x.g)G 0 (x,£. W ) , (3.9) 

where j3(x,£) is the angle between the normal vector i>(£) and v(x,£). 
Substituting (3.9) inside (3.1), and taking into account the definition (2.15) 
of F(x,k), this gives the following expression for the extrapolated field 

P e (x,<o) = kSfa) J ds F(f,k)n 0 (£)cosj3(x^)G 0 >t (x,£, W ) (3.10) 

which will be used in next section. 

It is worth noting that the definition (3.1) of the backpropagated field 
differs from the definition selected in [14] and [16], where P e (x,«) was given 
by the Kirchhoff integral 

P e (x, W ) = W(<o) J* ds[P s (£,«)v s G 0 ,t (x,f,(o) 

" v £ P s (f, U )G 0 X (x,g,(o)].u(£) . (3.11) 

The motivation for considering the expression (3.1) instead of (3.8) is that 
it is simpler, so that our subsequent derivations will be easier to follow. 
There is however no practical reason why (3.1) should be used instead of 
(3.11). An important difference between (3.1) and (3.11) is that the 
definition (3.1) requires only the knowledge of the scattered field P s (£, W ) 
along the receiver array T, whereas (3.11) requires also the knowledge of the 
normal derivative ^ P s (£,w) along T. Since this derivative cannot usually be 
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measured, it may appear at first sight that the definition (3.11) cannot be 

implemented directly. However, as was already observed above, it was shown by 

Esmersoy [14] that with the definition (3.11), P e (x, t) can be computed in the 

time domain by using the homogeneous equation (3.2) with the boundary 

condition (3.4), so that the normal derivative 3 P (£,«) is not needed. 

3 n s 

An argument which indicates that the definitions (3.1) and (3.11) are 

approximately equivalent is as follows. First, substitute the approximations 

r ik<f> (£) 

P s (£,m) = k S a s (£)e S (3.12a) 

G 0 (x,£,u) = k °a(x,£)e ik,p (-'£) (3.12b) 

inside (3.11), and take into account (3.9), so that 

r +1 

3 P (£,«) = ik S n 0 (£) cos0 s (£)P s (£,a,) (3.13a) 

dn 

r 0 +1 

§ G 0 (x,£.«) = ik n 0 (£) cosP(x,£)G 0 (x,£,w) (3.13b) 

3n 

where /3 (£) and /3(x,£) denote respectively the angles that v,-<p (£) and 
s £ s 

Vg<p(x,£) make with the normal to F. The integrals along T of the first and 
second term in (3.11) are now expressed as 

Ej = -ik r J ds nQ(£)a s (£)a*(x,|[)cos 0 (x,£)e ik ^s^ *(*,£)) (3.14a) 

E 2 = ” lkr J I dS n 0^^ a s^^ aX ^-’£^ cosP s^^ elk ^ S ^^ (3.14b) 

where r = Tq + r g + 1. The only difference between these two integrals is the 
weighting factors cosj3(x,£) and cos/3 g (£) which appear in (3.14a) and (3.14b) 
respectively. If the method of stationary phase is used to evaluate these 
integrals, we find that the stationary points are given by 

( V £<P S (£) “ Vg<Kx,£)).t(£) = 0 


(3.15) 
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where t(£) denotes the unit vector tangent to the array F at point f. But 
v^ s (f) and Vg<p(x,£) satisfy the eikonal equation and have therefore the same 

magnitude n^fg). In addition, they are both oriented in the outward direction 
with respect to F, so that at a stationary point we must have 

Y> s (£) = Y>(x,£) (3.16) 

which in turn implies that 

cos P s (l ) = cos/3(x,£) (3.17) 

Thus, at stationary points of (3.14a) and (3.14b) the weighting functions 

appearing in the two integrals are equal, so that Ej and Emake equal 

contributions to the leading order asymptotic expansion of P (x,w). 

e — ' 

Consequently, except for a factor 2 which can be incorporated in W(w), the 
expressions (3.1) and (3.11) are approximately the same. Note, however, that 
our analysis is only approximate since it relies on high frequency 
asymptotics. It turns out that there exist several special cases for which 
the two terms in (3.11) are exactly equal. This is the case for example when 
the background medium is constant and the receiver T is a straight line [19]. 
4. Inversion Method 
4.1 Reconstruction Technique 

The inversion procedure that we propose is just an extension of a method 
which was developed earlier by Esmersoy [14] for solving the linearized 
inversion problem for a point source with a constant background. The first 
step in our method is, for every point x in the medium, to image the 
backpropagated field P e (x,t) at the source travel time t = <p(x, tj)/c. Here 

represents the time needed for waves originating from the source 77 to 
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reach x. This gives the image 

°(x) = P e (?, | <f>(x,g)) = ^ f dw p e (x,w)e _1 c (4.1) 

The image 0(x) can be viewed as the migrated image obtained by applying the 
source travel-time imaging rule to the extrapolated field (see [26] for the 
description of a migration technique based on this principle). 

Then, changing the integration variable from w to k=w/c, taking into 
account the representation (3.10) for P e (x,«) and identities (2.13), (2.15), 
and using the geometrical optics approximation for Gq(*,*,«), we obtain 

0(x) = I dx’ N(x,x’)U(x’) , (4.2) 

J v 

where 

N(x.x') = J dk k 3/2 e 11 r/ 4 W(oi) J ds n Q (f )cosP(x.f) . 

• *(x\£) . ,4. 3) 

Let now 

“(S) - V £ °W*V(*’H> (4-4) 

be the reconstructed value of the potential U(x). Since the vector V x <p(x,i]) 

is tangent to the ray linking the source 3 to point x, and has magnitude 

n Q(x), this value is obtained by taking the gradient of the migrated image 

0 (x) in the direction of the ray linking 3 to x, and by scaling the resulting 

expression with n Q (x)/a(x, 3 ). The scaling by n Q (x)/a(x, 3 ) is needed to take 

dispersion effects into account. Then, from (4.2) we find that 

U(x) = T dx’M(x,x’)U(x’) (4.5) 

J V 

where 

=S(i^)' v x N( 5-5 , )*v s »(x.!!) • 


(4.6) 
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Our objective in the remainder of this section will be to obtain a simple 
representation for M(x.x’), which will be used to show that the reconstructed 
potential U(x) is a good approximation of U(x). 

4.2 Representation of M fx.x’l 

The starting point in our derivation of a representation for M(x,x’) is 
the expression (4.6), where N(x.x’) is given by (4.3). To evaluate v x N(x,x’), 

we assume that 

|v x (cos /3(x,g)a(x,g)) | « |kcos 0(x,£)a(x,£)v x <I>(x,£) | , (4.7) 

which, again, is an approximation of the same type as the geometrical optics 
approximation. This gives 

V x N(x,x').V x <p(x,r}) * - Ip dk k 5/2 e i 7 r/ 4 W(w) J ds n Q (£) cos P(x.£) 

• a ( 5’ >£) a (x,f)v x $(x,f)*v x <p(x,t])e lk ^^- ’£) . (4.8) 

Next, we note that 

= u(x,£) + u.(x) , (4.9) 


where the vectors 


A 

u 4 (x) = v x <f(x, 2 ) 


A 

u (x<D = v x «p(x.£) 


(4.10a) 

(4.10b) 


are tangent to the rays linking 5 to x and £ to x, respectively, as shown in 
Fig. 2. These vectors are such that 

|u.(x)| = |u(x,£)| = n Q (x) , (4.11) 

and in the following the angular arguments of u.(x) and u(x,£) will be denoted 
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by 0.(x) and 0(x,£), respectively. Thus 

v x $ (?S>£)* v x < # , (x*Il) = n o(x)(! + cosa(x,£)), (4.12) 

A 

where a(x,£) = 0(x,£) - 0.(x) is the angle between the vectors u(x,£) and 
u i (x). Then, by combining (4.6), (4.8) and (4.12), we obtain 

M(x,x*) = " ^(x~g) dk k 5/2 e i7r/4 W(«) J ds n Q (£) cos J3(x,£). 

(1 + cosa(x,£)) A(x',f)a(x,£)e lk[$ ^ ,, £ ) " ^.f)] . (4.13) 

This expression can be simplified further if, following Beylkin [1], we 
make the following approximation 

A(x’,£) ^ A(x,£) (4.14a) 

$(x‘,£) ^ $(x,£) + v x $(x,£)*(x'-x) (4.14b) 

for x, x' e V and £ e T. A rigorous analysis of this approximation in terms 
of pseudodifferential operators was proposed by Beylkin. In this context, it 
was shown that (4.14) has the effect of retaining the singular component of 
M(x.x’), which is nonzero only when x* is in the vicinity of x, and of 
neglecting the smooth components of M(x,x’). When the smooth components of 
M(£’£’) are dropped from identity (4.5), only the rapid space variations of 

A 

U(x), such as the location and size of discontinuities, can be related to 
those of U(x). Equivalently, in the Fourier domain, only the high wavenumber 
part of the Fourier transform of U(x) is related to that of U(x). Another 
interpretation of appromixation (4.14) which is perhaps more appealing from a 
physical point of view consists in observing that when V is in the far field 
of the receiver array f, the distance |x - x’| between points x, x’eV is small 
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with respect to |x - £| where geT, and in this case the Taylor series 
expansion (4.14) is justified. Note however that for geophysical surveys, it 
cannot always be assumed that V is in the field of F, and in this case the 
justification of (4.14) proposed by Beylkin is more appropriate. 

The end effect of (4.14) is that (4.13) takes the form 
cn 2 (x) 

M(x.x') =- 1~~ - J*° dk k 5/2 e iir/4 W(w) J ds n Q (g)cos P(x,f)* 

ikv $(x,£)*(x’-x) 

(1 + cos a(x,g)) a (x,£)e - . (4.15) 

Now, consider an infinitesimal ray tube originating from x and centered around 
the ray linking x to ( e F, as shown in Fig. 3. If £’is an arbitrary point 
along the ray linking x to (, and if do and do' are the cross-sections of the 
tube at £ and £' respectively, we have [27] 

n 0 (£) a2 (x>£) dcr = n 0 (D a2 (x.£’)dCT* • (4.16) 

But da can be expressed in terms of the element ds of arc length along F as 

da = ds cos /3(x,£) . (4.17) 

Furthermore, if d0 is the angular span of the tube at x, and if £' is located 
at a very small distance p from x, we have 

n 0 (D ~ n 0 (?S) (4.18a) 

da* = pd0 (4.18b) 

a2( S'S') ~ 8 ttt1q(x)p < 4 ' 1Sc > 

where to obtain (4.18c) we have used the asymptotic form (2.10) of a 
homogeneous Green’s function, and the fact that in the vicinity of x, the 
medium is locally homogeneous. This gives 

n 0 (£) a2 (5-£) cos P(x.£) ds = d0 , 


(4.19) 
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so that M(x,x') can be rewritten as 

2 
ci 

M( 


cn 2 (x) 

x,x') =. .- I dk k^ 2 e i7r/, ^W(w) f ds (1 + cos a(x,£)). 

16tt J I 


de , c . ikv «>Cs £)-(s -s) 

* d? - 

Now, consider the change of variables 

T(x) : (s,k) -» g = kv ^(x,f(s)) 


(4.20) 


(4.21) 


and let C(x) be the image of I x 18 under this transformation. Since g 
depends linearly on k, the domain C(x) is a cone whose span varies with x. It 
was shown by Beylkin [1] that the Jacobian of this transformation is given by 
J(x, s,k) = |k|n 2 (x)(l + cos a(x,£(s))) (x,£(s)) . (4.22) 

Consequently, if we select 


4 -i7r/4 

,M = - 7^ ’ (423) 

the weighting function M(x.x') can be expressed as 

M(x,x') = —^— 5 - f dg e*^ , (4.24) 

(2ttj J C(x) 

which is the desired representation for M(x.x'). 

From (4.24), we see that when M(x.x') is viewed as a function of d = x'-x 
parametrized by x, it is the inverse Fourier transform of the characteristic 
function 


(l for g e C(x) 

X (X,g) =! * (4.25) 

( 0 otherwise 

This shows that M(x,x') is a partial reconstruction of an impulse, where the 

cone C(x) specifies the range of available projections. In the special case 
2 

when C(x) = IR , which corresponds to an experiment geometry where the 
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receiver array T surrounds V, we have M(x,x') = 6(x-x'), so that in this case 
we can conclude that 

U(x) = U(x). (4.26) 

The identity (4.26) seems to suggest that our inversion procedure 
recovers exactly the scattering potential U(x). It is worth remembering at 
this point that identity (4.26), as well as the representation (4.24) for the 
kernel M(x,x’) are only valid within the limits imposed by the approximations 
we have made. These approximations are of course the Born approximation, but 
also the geometrical optics expansion (2.8) and approximation (4.14). It 
turns out that these last two approximations impose some very severe 
restrictions on our interpretation of (4.24) and (4.26). As we observed 
above, approximation (4.14) is automatically satisfied if V is in the far 
field of T. Otherwise, as was shown by Beylkin [1], the smooth components of 
M(x.x’) are neglected, so that expression (4.25) for the Fourier transform of 
M(x,x') is only valid when the wavevector g is large, and therefore identity 
(4.26) should be interpreted in the Fourier domain as 

U(E) = U(E) (4.27) 

for large g. 

It will be shown now that the geometrical optics expansion (2.8) imposes 
similar restrictions on the validity of relations (4.24) and (4.26). To see 
this, note that since k has been assumed large, there exists some constant r 
such that |k|> r. Then the transform variable g defined by (4.21) takes 
values over a set C r (x) which is the image of I x {k:|>r} under T(x). This 
set is obtained by subtracting from the cone C(x) a region in the vicinity of 
g = 0. The set C r (x) is characterized in detail in section 5. The main 
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aspect of this characterization is that although k is bounded away from zero, 
the magnitude of p = |g| of the wavevector g is not necessarily nonzero, since 
it depends on the length of the vector V x $(x,£) given by (4.9). In fact, when 

the receiver £ is located on the other side of V with respect to the source tj 
which may occur for a vertical seismic profiling experiment where the source 
is on the surface of the earth and the receivers along a vertical borehole, 
the length of V x $(x,£) is close to zero. Nevertheless, the constraint |k|>r 

has the effect that the Fourier transform relation (4.25) for M(x,x’) is 
restricted to g e C r (x). In the case of complete receiver coverage, this also 
implies that equality (4.27) between the Fourier transforms of U(x) and U(x) 
holds primarily for large values of g. 

Thus, both approximations (4.14) and the high frequency asymptotics that 
we have employed have the effect of restricting the validity of our inversion 
method to the high wavenumber region. 

To interpret the filter W(w) which is used to process the scattered field 
P (£.«). also observe that since k = w/c, we can rewrite 

W(co) = 4c 1/2 —- - X — T7 k . (4.28) 

(-i«) (-i«) 

so that W(u) corresponds to an integration followed by a square-root 
integration. 

4.3 Summary 

To conclude, let us review our reconstruction procedure step by step. 

1) First, the scattered field P s (£,w) observed at the receivers is 
filtered with W(w), which requires performing an integration, followed by a 


square-root integration. 
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2) The receivers are replaced by doublet sources, and the filtered time 
traces are used as source wavelets. Then, for this distribution of sources, 
the backpropagated field p e (x,t) is computed by a finite-difference scheme. 

3) The backpropagated field P e (x,t) is imaged at the source travel time 
<f>(x,ij)/c. This gives the migrated image 0(x). 

4) The reconstructed potential U(x) is given by (4.4), which is obtained 
by taking the gradient of 0(x) along the ray linking the source tj to point x, 
and by scaling the resulting expression with nQ(x)/a(x,£). 

The inversion procedure described above requires the computation of the 
extrapolated field p (x,t) and the evaluation of the phase <p(x,n) and 
amplitude a(x,i]) for all points in the medium. Since the location tj of the 
source is fixed, this scheme requires the use of a ray tracing algorithm only 
for the source rj. By comparison, the inverse GRT [l]-[3] requires the use of 
a ray tracing procedure not only for tj, but also for every receiver f located 
along F. Our method can therefore be viewed as having replaced the use of a 
ray tracing scheme for receivers along F by the computation of p e (x,t), which 
can be done in one batch operation, instead of receiver by receiver. The 
advantages and disadvantages of our inversion method with respect to the 
inverse GRT are primarily those of finite difference schemes with respect to 
ray tracing methods. Thus, when the number of receivers is large, and when 
the finite-difference scheme which is used to compute P e (x,t) does not require 
too many grid points, our reconstruction technique is likely to be faster than 
the inverse GRT. 

5. Characterization of the Cone C(x) 

The accuracy of the reconstruction method which has been obtained in the 
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previous section depends in a crucial way on the angular aperture of the cone 
C(x) appearing the representation (4.24) of the weighting function M(x.x'). 
Indeed, as the angular aperture of C(x) increases, within the approximations 
we have made M(x.x') gets closer to an impulse S(x-x'), and the Fourier 
transform U(g) of the reconstructed potential becomes a better approximation 
of the U(g) in the high wavenumber region. We will show now that the angular 
range of C(x) is purely a factor of the experiment geometry, and does not 
depend on the reconstruction technique that we have employed. In fact, as was 
already observed above, the same cone C(x) appears also in the inverse GRT. 
Specifically, it will be shown that the angular range of C(x) depends on the 
degree of coverage of the domain V which is provided by the receiver array T, 
and on the relative location of the source rj with respect to point x. 

The first step is to note from (4.21) that v x $(x,£) indicates the 

direction of a ray inside the cone C(x), so that as £ varies along T, V x $(x,£) 
spans the cone C(x), as shown in Fig. 4. The relation (4.9) also shows that 

A A 

v x $ (*'£) is the sum the two vectors u(x,£) and u^x), which have the same 
length, and have angular arguments 0(x,£) and 0.(x), respectively. Since 

A A 

u(x,£) and u^x) have the same length, the angle separating these two vectors 
is bisected by V x $(x,£), as indicated in Fig. 5, so that the angular argument 

^(x,£) of V x $(x,£) is given by 

*(x.£) = | (0(x,£) + 0 i (x)) . (5.1) 

Furthermore, since the source tj is fixed, u^x) and therefore O^x) are fixed, 
so that in (5.1) only 0(x,£) varies as £ moves along T. Consequently, if we 
assume that 0(x,£) varies over the angular range «(x) = [O^x), @ 2 (x)] as £ 
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moves along F, the angular span of the cone C(x) is 

#(X) = [| (e^x) + 0.(x)), | (0 2 (x) + 0.(x))] . (5.2) 

l 

An interesting feature of this result is that the total aperture of C(x) is 2 
(0 2 (x)-0^(x)), which is half the range of ®(x). 

To interpret the above result, consider the special case when the 
background medium is constant . In this case, it is possible to determine more 
precisely the effect of the receiver array F on the span #(x) of C(x). To do 
so, we will use for the array T a model of the type considered by Porter [20] 
and Esmersoy and Levy [16], where it is assumed that T is asymptotic to two 
lines wiht angles and a 2 with respect to the horizontal axis, as shown in 
Fig. 6. In exploration geophysics, this model covers the case when the 
scattered field is measured by receivers located on the surface of the earth 
(ctj =0, a 2 = tt) , or the offset vertical seismic profiling geometry, where the 
receivers are located on the surface of the earth and along a vertical 
borehole, say to the right of domain V (a^ = g-, a 2 = ir), or even the case when 
the receivers are located above V and along two vertical boreholes on both 
sides of V (a^ = a 2 = |r) • 

Then, since the background medium is constant, the rays linking a point x 
in the medium to the source ij and to receivers along T are straight lines . 
Thus, the angle 0^(x) is the angle of the line from tj to point x, and by 
keeping track of the lines linking points along T to x, we see that as £ moves 
along F, u(x,g) sweeps a domain D illustrated in Fig. 6, whose angular range 
is 0 = [ag-TT.a^+ir]. This angular range is independent of x and is completely 
parameterized by the angles a^ and a 2 describing the array F. Consequently, 
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the range 

*(*) = [| (a 2 -T+0 i (x)) , I (a 1 -Hr+0 i (x))] (5.3) 

of cone C(x) is purely a factor of the experiment geometry, and depends on x 
only through the angle 0,.(x) describing the relative location of the source g 
and point x. 

The last point that needs to be clarified is the effect of the high 
frequency constraint |k|>r on the range of values of the wavevector p under 
the transformation (4.21). As indicated at the end of section 4, in this case 
g & C r (x) where C r (x) is obtained by subtracting from the cone C(x) a set 
located in the vicinity of g = 0. To characterize this set, note from (4.21) 
that 

p = |g| = |k| |v x <?>(x,£) | , (5.4) 

and taking into account the fact that the vector V x $(x,£) is the sum of the 

/V A 

vectors u(x,£) and u^x) which have the same length n Q (x), we find that 

l v x $ (x,£)l = 2n 0 (x) cos((0(x,f)-0 i (x))/2) (5.5) 

Thus, in the direction ^(x,g) = (0(x,g)+0 1 (x))/2, the wavevector g must be 
such that 

P > 2rn 0 (x) cos((0(x,£)-0.(x))/2) = 2rn Q (x) cos(^0.(x)) . (5.6) 

But p = 2rn^(x) cos^-O^x)) is the equation of a circle centered at rUjfx) 
and of radius rn Q (x). This shows that C r (x) is obtained by subtracting from 
the cone C(x) the points which are inside a disk of radius rn Q (x) centered at 
ru -[(x), or which are inside the symmetric image of this disk with respect to 
the origin, as indicated in Fig. 7. 

An interesting aspect of the above result is that when ^ = 0.(x)+ir, which 

1 2 
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corresponds to the case where 0(x,£) = 0.(x)+tt, or equivalently 
u(x,£) = -u.(x), then the origin 0 belongs to the domain C r (x). This 
indicates that some coverage in the low wavenumber region is possible for 
certain source-receiver geometries. Note that when u(x,£) = -u..(x), the 
receiver £ is located on the other side of V with respect to the source 77 , but 
on the same ray. This geometry is of a tomographic nature, and arises in 
exploration geophysics for vertical seismic profiling or borehole to borehole 
surveys, where in the last case the source and receivers are located in two 
different vertical boreholes. By comparison, note that when 'p = 0^(x), which 

A A 

corresponds to 0(x,£) = 0,.(x) and u(x,£) = -u^(x), the wavenumber p must be 
such that p > 2rnQ(x), and g is therefore bounded away from zero. This 
special case corresponds to the so-called zero-offset experiment geometry, 
where the source and receiver coincide. Such a geometry can be used as a 
model for surface surveys in exploration geophysics, where the source and 
receivers are located on the surface of the earth at relatively small distance 
from each other. The above observations indicate that in order to gain some 
information about the slow space variations of the scattering potential U(x), 
a tomographic experiment geometry must be employed. 
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6. Conclusion 

In this paper, we have obtained a new solution to the variable background 
Born inversion problem for the case when the medium is probed by a single 
point source. This solution generalizes a reconstruction technique proposed 
earlier by Esmersoy [14] for the constant background case. In our 
reconstruction method, the backpropagated field p e (x,t) is first computed by a 
finite difference sheme, and is imaged at the source travel times. This gives 
the migrated image 0(x), and by taking the gradient of this image along rays 
linking the source to every point in the medium, and scaling the resulting 
image appropriately, we obtain the reconstructed potential U(x). When the 
receiver array provides a total coverage of the domain where the scattering 
potential U(x) is concentrated, the reconstructed potential recovers exactly 
the rapid space variations of U(x), or equivalently the high wavenumber part 
of its Fourier transform. In general, within several approximations which 
restrict the validity of our reconstruction procedure to the high wavenumber 
region, it is shown that U(x) corresponds to the convolution at point x of 
U(*) with a weighting function obtained by reconstructing an impulse from its 
projections over a cone C(x). The angular range of the cone C(x) depends only 
on the experiment geometry, and not on the reconstruction method that we have 
employed. 

The most significant computational requirement of the procedure described 
above is the computation of the extrapolated field p e (x,t) with a 
finite-difference scheme. By comparison, the generalized Radon transform 
inversion method [l]-[3] requires the use of a ray tracing scheme, and the 
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evaluation of the phase and amplitude along each ray, for every receiver along 
the receiver array T. We expect therefore that, for a wide class of problems, 
our method will be faster than the inverse GRT but a detailed comparison of 
the relative efficiency of these two techniques will be necessary. A number 
of factors are likely to influence this comparison. The first factor is the 
number of receivers along T- clearly, as this number grows, the numerical 
complexity of the inverse GRT increases significantly, whereas the number of 
operations required to compute the backpropagated field P e (x,t) remains 
approximately the same. Another factor is the complexity of the background 
refraction index profile n^f*). It is known that ray tracing schemes perform 
very well for simple media, but are relatively inaccurate, and difficult to 
use for complex media. By comparison, the finite-difference method is a brute 
force technique which is not affected significantly by the complexity of the 
background profile n^f*). Finally, another important factor in our comparison 
will be the relative location of the domain V that we want to image with 
respect to the receiver array T. If this domain is not too far away from T, 
the finite-difference technique is relatively easy to use, since it does not 
require the computation of p e (x,t) over a very large region of space. 

However, if V is far away from F, the finite-difference technique becomes very 
slow, whereas the ray tracing approach is not significantly affected by the 
increase in size of the domain that we consider. The above observations seem 
to suggest that the inverse GRT and the method that we propose here are almost 
complementary, in the sense that one method will perform best on problems for 
which the other is least suited. However, these conclusions are rather 
tentative, since as mentioned above, a detailed comparison of these two 
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methods has yet to be performed. Note that our reconstruction technique gives 
good results in the constant background case [14]. 

Several assumptions have been made throughout this paper. The first of 
these is that the medium that we consider is two-dimensional. The 
reconstruction procedure which has been presented above can be extended in a 
straightforward way to 3-D media. Another significant restriction which has 
been imposed is that the background refraction index n^(*) is a smooth 
function. This assumption is somewhat unrealistic in practice, since most 
geological structures exhibit significant discontinuities. It would therefore 
be of interest to extend our reconstruction technique to the case when the 
background index n^f*) is discontinuous. The main difficulty in this context 
is that at interfaces where n^f*) is discontinuous, the incident waves are 
partly reflected, so that multiply reflected waves exist, and in this case the 
asymptotic form (2.8) of the Green’s function is not valid. If the reflected 
waves are neglected, and if only the transmission losses at interfaces are 
taken into account, it was shown in [1] how to reconstruct U(•). However, it 
would be ultimately desirable to obtain a reconstruction technique that takes 
into account the reflected waves appearing in the background model. In 
addition, we note that our reconstruction technique applies only to before 
stack data, i.e. to data collected from a single experiment. In actual 
geophysical surveys, several experiments are carried out for different source 
and receiver locations. It would therefore be useful to find a way of 
combining the reconstructed potentials Ikfx), 1 < i < k obtained for several 
experiments, where k is the number of experiments. Note that in this case the 
cones C^(x), 1 < i < k at a point x will provide different angular coverages 
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for different experiment geometries. A simple averaging of the potentials 
U.(x) does not give a satisfactory solution to this problem, since this 
average would have for effect to weight too heavily the regions where the 
cones (L(x) overlap. Another way of formulating this problem is to consider 
the zero-offset geometry, where the data is given for coincident 
source-receiver pairs located along a curved array. This data is obtained by 
applying the common depth point (CDP) stacking process [23] to the data 
collected from a large number of experiments. For a 2 1/2-D zero-offset 
geometry, and when the background refraction index is constant, an inversion 
technique which relies on the backpropagated field was proposed by Esmersoy 
[14]. We expect that this inversion technique can be extended to the variable 
background case by using ideas similar to those that have been discussed in 
this paper. 
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Figure Captions 

Figure 1: Scattering experiment. The medium is probed by an impulsive point 
source located at tj and the scattered field is observed along the receiver 
array T. 

Figure 2- Parametrization of the rays linking point x to the source tj and to 
receiver £. The doublet sources generating the backpropated field are 
indicated by + and - signs. 

Figure 3: Infinitesimal ray tube originating from point x and centered around 
the ray linking x to receiver £. 

Figure 4: Cone C(x) spanned by V x $(x,£) as £ moves along T. 

Figure 5: Construction of V x $(x,£) by summation of the vectors u(x,f) and 
u.(x). 

Figure 6: Receiver array with singular aperture [a^, a^]. D denotes the 

A 

domain swept by the vector u(x,£) as £ moves along T. 

Figure 7: Range C r (x) of the wavevector g under the high frequency constraint 

|k|>r. 
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